/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TAMOMDIS_M
  \brief Calculate core momentum distribution with angular momentum component m
  specified. This is a class to assist class TAMOMDIS.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/09
  \date Last modified: 2022/04/05 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <cstring>
#include <iostream>
#include <iomanip>
#include <fstream>
#include "TAMOMDIS_M.h"
#include "TAMOMDIS.h"
#include "TAFun.h"
#include "TAMath.h"
#include "TABound.h"
#include "TABessel.h"
#include "TASMatrix.h"
#include "TAIntegrate.h"
#include "TALegendre.h"

using std::cout;
using std::endl;
using std::flush;
using std::norm;
using std::conj;
using std::ofstream;
using std::setw;

#define sqr TAMath::sqr
#define clm2 TAMath::Clm2
#define romb TAIntegrate<>::Romberg
#define rombc TAIntegrate<cdouble>::Romberg
#define simp TAIntegrate<>::Simpson
#define simpc TAIntegrate<cdouble>::Simpson
#define leg TALegendre::Legendre

static auto fn = TAMath::factln; // fn(n) = ln(n!)
static constexpr double HBARC = TAMath::hbarc();
static constexpr double PI = TAMath::Pi(), TWOPI = 2.*PI, FOURPI = 4.*PI;
static constexpr cdouble I(0., 1.); // the imaginary number unit

// all in fm
const double TAMOMDIS_M::RMAX = 30.;
const double TAMOMDIS_M::BMAX = 30.;
const double TAMOMDIS_M::ZERO = 1.e-5;
const double TAMOMDIS_M::MOMMAX = 500./HBARC; // the maximum momentum

TAMOMDIS_M::TAMOMDIS_M(TAMOMDIS *mom) : fMOM(mom){
  fSc = fMOM->Sc();
  fSv = fMOM->Sv();
  fRL = fSc->GetRL();
  if(BMAX*TAMath::Sqrt2() > fSc->BMax() || BMAX*TAMath::Sqrt2() > fSv->BMax())
    TAException::Error("TAMOMDIS_M", "ctor: BMAX too large, cannot be larger than s-matrice's BMAX");
  l = fMOM->GetBound()->Getl();
  if(l < 0 || l > 10) TAException::Warn("TAMOMDIS_M", "ctor: abnormal l: %d", l);

  const int ng = fNG = fSc->GetNGaus(); // number of gaussians used in the sum-gaus fit
  fAlphajR = new double[ng]{}; fAlphajI = new double[ng]{};
  // memcpy(fAlphajR, fSc->GetAlphajR(), sizeof(double)*ng); // complete the gaussian ->
  // memcpy(fAlphajI, fSc->GetAlphajI(), sizeof(double)*ng); // -> expansion fit
  fBetaM2 = new double[ng];
  for(int j = ng; j--;) fBetaM2[j] = sqr(j/fRL); // 1./(RL/j)^2
  // assign array fMomArr; the HBARC is for transform from MeV/c -> fm^-1
  const double h = MOMMAX/(kNmom-1); // MeV/c per nucleon -> fm^-1
  fMomArr = new double[kNmom];
  for(int i = kNmom; i--;) fMomArr[i] = h*i; // element with i = 0 stores the total c.s.
  // set bv and rho mesh //
  dr = (RMAX-ZERO) / (NR-1); dr2 = RMAX*2. / (NR2-1);
  db = (BMAX-ZERO) / (NB-1);
  dphi = PI / (NPHI-1); dphi2 = TWOPI / (NPHI2-1);
  for(int i = 0; i < NB; i++) bv[i] = db*i + ZERO;
  for(int i = 0; i < NR; i++) r[i] = dr*i + ZERO;
  for(int i = 0; i < NR2; i++) r2[i] = dr2*i - RMAX + ZERO;
  for(int i = 0; i < NPHI; i++) phi[i] = dphi*i;
  for(int i = 0; i < NPHI2; i++) phi2[i] = dphi2*i;
} // end the constructor

TAMOMDIS_M::~TAMOMDIS_M(){
  if(fAlphajR) delete [] fAlphajR;
  if(fAlphajI) delete [] fAlphajI;
  if(fBetaM2) delete [] fBetaM2;
  if(fMomArr) delete [] fMomArr;
} // end the destructor

double TAMOMDIS_M::ParallelStr(int m, double *momStr){ // accumulate the results
  m = abs(m);
  if(m > l) TAException::Error("TAMOMDIS_M", "ParallelStr: illegal m=%d", m);

  cout << "\033[33;1mCalculating stripping mom-dis for m = " << m << "...\033[0m" << endl;
  cout << "Integrating Sc over phi. Please wait. This may take a while..." << endl;
  // assign blm2 //
  for(int i = 0; i < NB; i++){
    for(int j = 0; j < NR; j++){
      if(!fBLMSQ[i][j]){ // assign fBLMSQ only once for all m
        const auto &intphi = [=](double phi){ // integrate |sc|^2 over phi
          return Sc2(sqrt(fabs(bv[i]*bv[i]+r[j]*r[j]-2.*bv[i]*r[j]*cos(phi))));
        }; // end lambda intphi
        fBLMSQ[i][j] = romb(intphi, 0., PI)*2.; // *2.: using symmetry at phi=pi
        // blm2[i][j] = Blm2(m, bv[i], r[j])*TWOPI; // note that this statement is problematic
      } // end if
    } // end for over j
    cout << "b: " << setw(13) << bv[i] << " -> b = " << BMAX << '\r' << flush;
  } // end for over i

  cout << "Calculating stripping mom-distri" << endl;
  for(int i = 0; i < kNmom; i++) momStr[i] = ParallelStr(m, fMomArr[i]);
  // output to a text file //
  char file[256] = "";
  sprintf(file, "vis/strm%d.txt", m);
  ofstream fout(file);
  for(int i = kNmom; i--;) fout << -fMomArr[i]*HBARC << " " << momStr[i] << endl;
  for(int i = 1; i < kNmom; i++) fout << fMomArr[i]*HBARC << " " << momStr[i] << endl;
  fout.close();
  cout << "Stripping distribution has been written to file " << file << endl;
  return simp(kNmom, fMomArr, momStr)*2.*HBARC; // *2.: symmetry at kz=0
} // end member function ParallelStr

/// the integral sum_p{|B_lmp|^2} with ZINTEG. Ref[PRC.70.034609-12, the appendix]
double TAMOMDIS_M::Blm2(int m, double bv, double rho) const{
  if(!fAlphajR[1] && !fAlphajI[1])
    TAException::Error("TAMOMDIS_M", "Blm2: Gaus Fit not implemented yet~");

  static const int PMAX = 5;
  cdouble blm = 0.;
  const double tbr = 2.*bv*rho, br2 = bv*bv+rho*rho, bmr2 = sqr(bv-rho);
  for(int j = 0; j < fNG; j++){
    const double x = tbr*fBetaM2[j];
    double c = 0.; // the sum of besselI terms
    if(x < 350.){ // (p ? 1. : 2.) is because that I_{-v} = I_v
      for(int p = -PMAX; p <= PMAX; p++) c += TABessel::BesselI(p, x)*(p ? 2. : 1.);
      c *= exp(-br2*fBetaM2[j]);
    } // end if
    else{ // overflow for besselI, using BsiEmx=I_p(x)*exp(-x) instead
      for(int p = -PMAX; p <= PMAX; p++) c += TABessel::BsiEmx(p, x)*(p ? 2. : 1.);
      c *= exp(-bmr2*fBetaM2[j]);
    } // end else
    blm += (fAlphajR[j]+I*fAlphajI[j])*c;
  } // end for over j
  return norm(blm);
} // end member function Blm2

// calculate stripping (inelastic) knockout cross sections for momentum kz: dsig/dkz
double TAMOMDIS_M::ParallelStr(int m, double kz) const{
  double zint2[NR]{}, zzint[NR]{}; // integral over z
  // do integral along the z-direction and store in zint2
  for(int i = 0; i < NR; i++){
    for(int j = 0; j < NR; j++){
      const double z = r[j];
      const double rr = sqrt(r[i]*r[i]+z*z);
      zzint[j] = Rl(rr)*leg(l, m, z/rr)*
        ((l+m) & 1 ? -sin(kz*z) : cos(kz*z)); // exp(-ikz*z);
    } // end for over j
    zint2[i] = sqr(2.*simp(NR, r, zzint)); // *2.: using symmetry around z=0
  } // end for over i
  double bint[NB]{}, blm2[NB][NR]{}; // integral over bv
  for(int i = 0; i < NB; i++){
    for(int j = 0; j < NR; j++) blm2[i][j] = r[j]*zint2[j]*fBLMSQ[i][j];
    bint[i] = bv[i]*(1.-Sv2(bv[i]))*simp(NR, r, blm2[i]);
  } // end for over i

  // exp(fn(l-m)-fn(l+m))/FOURPI = clm^2/(2l+1); *10: fm^2->mb
  // /HBARC is necessary if momstr is against p instead of k
  const double momstr = simp(NB, bv, bint)*exp(fn(l-m)-fn(l+m))/FOURPI*10./HBARC * (m ? 2 : 1);
  cout << "kz: " << setw(8) << kz*HBARC << " -> " << MOMMAX*HBARC;
  cout << ", momstr: " << momstr << '\r' << flush;
  return momstr;
} // end member function ParallelStr

/// the same as ParallelStr, but for diffraction dissociation
double TAMOMDIS_M::ParallelDiff(int m, double *momDiff){
  m = abs(m);
  if(m > l) TAException::Error("TAMOMDIS_M", "ParallelDiff: illegal m=%d", m);

  cout << "\033[33;1mCalculating diffractive mom-dis for m = " << m << "...\033[0m" << endl;
  cout << "Calculating SPINT, please wait..." << endl;
  // calculate zint2 //
  cout << "Doing the z-integral..." << endl;
  double zint[NR]{}, zzint[NR]{}, rzint0[NR]{};
  for(int i = 0; i < NR; i++){
    for(int j = 0; j < NR; j++){
      const double z = r[j];
      const double rr = sqrt(r[i]*r[i]+z*z);
      zzint[j] = sqr(Rl(rr)*leg(l, m, z/rr));
    } // end for over j
    zint[i] = 2.*simp(NR, r, zzint); // *2.: using symmetry around z=0
    rzint0[i] = r[i]*zint[i];
  } // end for over i
  // renormalize the wavefunction for more accurate c.s.
  const double cc = 1./(TWOPI*simp(NR, r, rzint0));
  // calcualte ilm
  const double mc = fSc->PM(), mv = fSv->PM(), mp = mc+mv;
  const double ratioc = mv/mp, ratiov = mc/mp;
  cdouble ilm[NB]{}, rzint[NR]{};
  // calcualte ilm as a function of b //
  cout << "Calculating ilm. This may take a while..." << endl;
  for(int i = 0; i < NB; i++){
    const double bb = bv[i], b2 = bb*bb; // bv here is actually b (in c.m. frame)
    for(int j = 0; j < NR; j++){
      if(0. == fSSINT[i][j]){ // assign fSSINT only once for all m
        const double rc = ratioc*r[j], rv = ratiov*r[j];
        const double rc2 = rc*rc, rv2 = rv*rv, brc2 = b2+rc2, brv2 = b2+rv2;
        const auto &intphi = [=](double phi){
          const double cphi = 2.*bb*cos(phi);
          const cdouble &sc = Sc(sqrt(fabs(brc2-rc*cphi)));
          const cdouble &sv = Sv(sqrt(fabs(brv2+rv*cphi)));
          return sc*sv;
        }; // end lambda intphi
        fSSINT[i][j] = rombc(intphi, 0., PI)*2.;
      } // end if
      rzint[j] = r[j]*zint[j]*conj(fSSINT[i][j]);
    } // end for over j
    ilm[i] = simpc(NR, r, rzint)*cc; // renormalize a little bit
    cout << "b: " << setw(13) << bv[i] << " -> b = 30" << '\r' << flush;
  } // end for over i
  // calculate fSPINT //
  cout << "Calculating fSPINT. Also time-consuming tasks..." << endl;
  for(int i = 0; i < NB; i++){
    const double bb = bv[i], b2 = bb*bb; // bv here is actually b (in c.m. frame)
    for(int j = 0; j < NR; j++){
      if(!fSSINT2[i][j]){
        const double rc = ratioc*r[j], rv = ratiov*r[j];
        const double rc2 = rc*rc, rv2 = rv*rv, brc2 = b2+rc2, brv2 = b2+rv2;
        const auto &intphi = [=](double phi){
          const double cphi = 2.*bb*cos(phi);
          const cdouble &sc = Sc(sqrt(fabs(brc2-rc*cphi)));
          const cdouble &sv = Sv(sqrt(fabs(brv2+rv*cphi)));
          return norm(sc*sv);
        }; // end lambda intphi
        fSSINT2[i][j] = romb(intphi, 0., PI)*2.;
      } // end if
      fSPINT[i][j] = (fSSINT2[i][j] - ilm[i]*fSSINT[i][j])*cc; // cc: renormalize a little bit
    } // end for over j
    cout << "b: " << setw(13) << bv[i] << " -> b = " << BMAX << '\r' << flush;
  } // end for over i

  cout << "Here the mom-diff calculation begins..." << endl;
  for(int i = 0; i < kNmom; i++) momDiff[i] = ParallelDiff(m, fMomArr[i]);
  // output to a text file //
  char file[256] = "";
  sprintf(file, "vis/diffm%d.txt", m);
  ofstream fout(file);
  for(int i = kNmom; i--;) fout << -fMomArr[i]*HBARC << " " << momDiff[i] << endl;
  for(int i = 1; i < kNmom; i++) fout << fMomArr[i]*HBARC << " " << momDiff[i] << endl;
  fout.close();
  cout << "Diffractive distribution has been written to file " << file << endl;
  return simp(kNmom, fMomArr, momDiff)*2.*HBARC; // symmetry at kz=0
} // end member function ParallelDiff

// calculate stripping (inelastic) knockout cross sections for momentum kz: dsig_diff/dkz
double TAMOMDIS_M::ParallelDiff(int m, double kz) const{
  double zint2[NR]{}, zzint[NR]{}; // integral over z
  // do integral along the z-direction and store in zint2
  for(int i = 0; i < NR; i++){
    for(int j = 0; j < NR; j++){
      const double z = r[j];
      const double rr = sqrt(r[i]*r[i]+z*z);
      zzint[j] = Rl(rr)*leg(l, m, z/rr)*
        ((l+m) & 1 ? -sin(kz*z) : cos(kz*z)); // exp(-ikz*z);
    } // end for over j
    zint2[i] = sqr(2.*simp(NR, r, zzint)); // *2.: using symmetry around z=0
  } // end for over i
  cdouble bint[NB]{};
  cdouble rzint[NR]{};
  for(int i = 0; i < NB; i++){
    for(int j = 0; j < NR; j++) rzint[j] = r[j]*zint2[j]*fSPINT[i][j];
    bint[i] = bv[i]*simpc(NR, r, rzint);
  } // end for over i
  // Clm^2/(2l+1); *10: fm^2->mb
  const cdouble momdiff = simpc(NB, bv, bint)/(2.*l+1.)*10./HBARC * double(m ? 2 : 1);
  cout << "kz: " << setw(8) << kz*HBARC << " -> " << MOMMAX*HBARC;
  cout << ",  momdiff: " << momdiff.real() << '\r' << flush;
  return momdiff.real();
} // end member function ParallelDiff

// total cross section for the stripping channel
double TAMOMDIS_M::CSStr(int m) const{
  m = abs(m);
  if(m > l) TAException::Error("TAMOMDIS_M", "CSStr: illegal m=%d", m);

  char name[] = "vis/str.txt";
  ofstream fout(name);
  cout << endl << "\033[32;1mCalculating  stripping  c.s. for m = ";
  cout << m << "...\033[0m" << endl;

  // set rho and b mesh //
  double zint[NR]{}, zzint[NR]{}, rzint[NR]{}, bint[NB]{}; // integral over z, phi and bv
  // do integral along the z-direction and store in fint
  for(int i = 0; i < NR; i++){
    for(int j = 0; j < NR; j++){
      const double z = r[j];
      const double rr = sqrt(r[i]*r[i]+z*z);
      zzint[j] = sqr(Rl(rr)*leg(l, m, z/rr));
      // cout << "z: " << z << " rho: " << r[i] << endl;
      // cout << "Rl(" << rr << "): " << Rl(rr) << endl;
      // cout << "leg(" << l << ", " << m << ", " << z/rr << "): " << leg(l, m, z/rr) << endl;
      // cout << "zzint: " << sqr(Rl(rr)*leg(l, m, z/rr)) << endl;
      // getchar();
    } // end for over j
    zint[i] = 2.*simp(NR, r, zzint); // *2.: using symmetry around z=0
  } // end for over i
  for(int i = 1; i < NB; i++){
    for(int j = 0; j < NR; j++){
      const auto &intphi = [=](double phi){ // integrate |sc|^2 over phi
        return Sc2(sqrt(fabs(bv[i]*bv[i]+r[j]*r[j]-2.*bv[i]*r[j]*cos(phi))));
      }; // end lambda intphi
      rzint[j] = r[j]*zint[j]*romb(intphi, 0., PI, 1.e-3)*2.; // using symmetry at phi=pi
      // cout << std::scientific;
      // cout << "bv=" << bv[i] << " rho=" << r[j] << endl;
      // cout << "zint=" << zint[j] << " sint=" << romb(intphi, 0., PI, 1.e-3)*2. << endl;
      // getchar(); // DEBUG
    } // end for over j
    bint[i] = bv[i]*(1.-Sv2(bv[i]))*simp(NR, r, rzint);
    fout << bv[i] << " " << bint[i] << endl;
    cout << "b: " << setw(8) << bv[i] << " bint: ";
    cout << setw(13) << bint[i];
    // cout << " 1.-Sv2(bv[i]): " << 1.-Sv2(bv[i]);
    // cout << endl; // DEBUG
    cout << '\r' << flush;
    // getchar(); // DEBUG
  } // end for over i
  fout.close();
  cout << "The bv-integral v.s. bv has been stored in " << name << endl;
  cout << " and can be inspected via `cd vis/ && root -l str.C`" << endl;
  cout << "(provided ROOT is installed)." << endl;

  return simp(NB, bv, bint)*exp(fn(l-m)-fn(l+m))*5. * (m ? 2 : 1); // *5 = *10/2
} // end member function CSStr

// #define EX_ // calculate m',m - exchange term or not in sigma_diff

// total cross section for the diffraction dissociation channel
double TAMOMDIS_M::CSDiff(int m) const{
  m = abs(m);
  if(m > l) TAException::Error("TAMOMDIS_M", "CSDiff: illegal m=%d", m);

#ifdef EX_
  // adding exchange-terms //
  double reex = 0.;
  for(int mm = -l; mm <= l; mm++){
    // mm+m: even for nonzero results
    if((mm+m) & 1) continue;
    if(mm != m) reex += CSDiffEx(m, mm);
    if(m && (mm != -m)) reex += CSDiffEx(-m, mm);
  } // end for
#endif

  const char name1[] = "vis/diff.txt", name2[] = "vis/zint.txt";
  ofstream fout(name1), foutz(name2);
  cout << endl << "\033[32;1mCalculating diffractive c.s. for m = ";
  cout << m << "...\033[0m" << endl;

  // set rho, b and phi mesh //
  double zint[NR]{}, zzint[NR]{}, bint[NB]{}; // integral over z and bv
  double rzint1[NR]{}, phint1[NPHI]{}; // integral of |1-Sc*Sv|^2
  cdouble rzint2[NR]{}, phint2[NPHI]{}; // integral of (1-Sc*Sv)
  double rzint0[NR]{}; // for renormalization of the wavefunction
  // // the left-out exchange terms - 2023.1.2 //
  // double zintex[NR]{}, rzintex[NR]{}, phintex[NR]{};
  // do integral along the z-direction and store in fint
  for(int i = 0; i < NR; i++){
    for(int j = 0; j < NR; j++){
      const double z = r[j];
      const double rr = sqrt(r[i]*r[i]+z*z);
      zzint[j] = sqr(Rl(rr)*leg(l, m, z/rr));
    } // end for over j
    zint[i] = 2.*simp(NR, r, zzint); // *2.: using symmetry around z=0
    rzint0[i] = r[i]*zint[i];
    foutz << r[i] << " " << zint[i] << endl;
  } // end for over i
  cout << "The z-integral v.s. rho has been stored in " << name2 << endl;
  foutz.close();
  // renormalize the wavefunction for more accurate c.s.
  const double cc = 1./(TWOPI*simp(NR, r, rzint0));
  const double mc = fSc->PM(), mv = fSv->PM(), mp = mc+mv;
  const double ratioc = mv/mp, ratiov = mc/mp;
  for(int i = 0; i < NB; i++){
    const double bb = bv[i], b2 = bb*bb; // bv here is actually b (in c.m. frame)
    for(int j = 0; j < NR; j++){
      const double rc = ratioc*r[j], rv = ratiov*r[j];
      const double rc2 = rc*rc, rv2 = rv*rv;
      for(int k = 0; k < NPHI; k++){
        const double cphi = 2.*bb*cos(phi[k]);
        const cdouble &sc = Sc(sqrt(fabs(b2+rc2-rc*cphi)));
        const cdouble &sv = Sv(sqrt(fabs(b2+rv2+rv*cphi)));
        phint2[k] = sc*sv;
        phint1[k] = norm(phint2[k]);
      } // end for over k
      // *2.: using symmetry around phi=pi
      rzint1[j] = r[j]*zint[j]*2.*simp(NPHI, phi, phint1);
      rzint2[j] = r[j]*zint[j]*2.*simpc(NPHI, phi, phint2);
    } // end for over j
    bint[i] = bv[i]*(simp(NR, r, rzint1) - norm(simpc(NR, r, rzint2))*cc);
    fout << bv[i] << " " << bint[i] << endl;
    cout << "b: " << setw(8) << bv[i] << " bint: ";
    cout << setw(13) << bint[i] << '\r' << flush;
  } // end for over i
  fout.close();
  cout << "The b-integral v.s. b has been stored in " << name1 << endl;
  cout << "They can be inspected via `cd vis/ && root -l str.C`" << endl;
  cout << "(provided ROOT is installed)." << endl;

  const double re = simp(NB, bv, bint)*TWOPI/(2.*l+1.)*cc*10. * (m ? 2 : 1); // *10: fm^2->mb

#ifdef EX_
  cout << "re: " << re << " mb, reex: " << reex << " mb" << endl; // DEBUG
  cout << "reex/re = " << reex/re << endl; // DEBUG
  // getchar(); // DEBUG
  return re + reex;
#else
  return re;
#endif
} // end member function CSDiff

// the left-out exchange term - 2023.1.2
double TAMOMDIS_M::CSDiffEx(int m, int mm) const{
  if(fabs(m) > l) TAException::Error("TAMOMDIS_M", "CSDiffEx: illegal m=%d", m);
  if(fabs(mm) > l) TAException::Error("TAMOMDIS_M", "CSDiffEx: illegal mp=%d", mm);
  if(mm == m) TAException::Error("TAMOMDIS_M", "CSDiffEx: mp==m occurred");
  const double dm = mm - m;

  const char name1[] = "vis/diffex.txt", name2[] = "vis/zintex.txt";
  ofstream fout(name1), foutz(name2);
  cout << endl << "\033[32;1mCalculating diffractive-ex c.s. for m = ";
  cout << m << " and mm = " << mm << "...\033[0m" << endl;

  // set rho, b and phi mesh //
  double zint[NR]{}, zzint[NR2]{}, bint[NB]{}; // integral over z and bv
  cdouble rzint2[NR]{}, phint2[NPHI]{}; // integral of (1-Sc*Sv)
  // do integral along the z-direction and store in fint
  for(int i = 0; i < NR; i++){
    for(int j = 0; j < NR2; j++){
      const double z = r2[j];
      const double rr = sqrt(r[i]*r[i]+z*z);
      const double ct = z/rr;
      zzint[j] = sqr(Rl(rr))*leg(l, m, ct)*leg(l, mm, ct);
    } // end for over j
    zint[i] = simp(NR2, r2, zzint);
    foutz << r[i] << " " << zint[i] << endl;
  } // end for over i
  cout << "The z-integral v.s. rho has been stored in " << name2 << endl;
  foutz.close();
  // renormalize the wavefunction for more accurate c.s.
  const double mc = fSc->PM(), mv = fSv->PM(), mp = mc+mv;
  const double ratioc = mv/mp, ratiov = mc/mp;
  for(int i = 0; i < NB; i++){
    const double bb = bv[i], b2 = bb*bb; // bv here is actually b (in c.m. frame)
    for(int j = 0; j < NR; j++){
      const double rc = ratioc*r[j], rv = ratiov*r[j];
      const double rc2 = rc*rc, rv2 = rv*rv;
      for(int k = 0; k < NPHI2; k++){
        const double cphi = 2.*bb*cos(phi2[k]);
        const cdouble &sc = Sc(sqrt(fabs(b2+rc2-rc*cphi)));
        const cdouble &sv = Sv(sqrt(fabs(b2+rv2+rv*cphi)));
        phint2[k] = (sc*sv)*exp(I*dm*phi2[k]);
        // phint2[k] = exp(I*dm*phi[k]);
      } // end for over k
      // *2.: using symmetry around phi=pi
      rzint2[j] = r[j]*zint[j]*simpc(NPHI2, phi2, phint2);
      // rzint2[j] = r[j]*zint[j]*2*simpc(NPHI, phi, phint2);
    } // end for over j
    bint[i] = bv[i]*(norm(simpc(NR, r, rzint2)));
    fout << bv[i] << " " << bint[i] << endl;
    cout << "b: " << setw(8) << bv[i] << " bint: ";
    cout << setw(13) << bint[i] << '\r' << flush;
  } // end for over i
  fout.close();
  cout << "The b-integral v.s. b has been stored in " << name1 << endl;
  cout << "They can be inspected via `cd vis/ && root -l str.C`" << endl;
  cout << "(provided ROOT is installed)." << endl;

  const double re = simp(NB, bv, bint)*TWOPI/(2.*l+1.)*10.*clm2(l,m)*clm2(l,mm);
  cout << "m = " << m << " mm = " << mm; // DEBUG
  cout << " CSDiffEx = " << re << " mb." << endl; // DEBUG
  // getchar(); // DEBUG
  return re; // *10: fm^2->mb
  // return simp(NB, bv, bint)*TWOPI/(2.*l+1.)*cc*10. * (m ? 2 : 1); // *10: fm^2->mb
} // end member function CSDiffEx


/// the radial wavefunction
double TAMOMDIS_M::Rl(double r) const{
  return fMOM->Rl(r);
} // end member function Rl

cdouble TAMOMDIS_M::Sc(double b) const{ return fSc->SMatrix(b); }
double TAMOMDIS_M::Scr(double b) const{ return fSc->SMatrixR(b); }
double TAMOMDIS_M::Sci(double b) const{ return fSc->SMatrixI(b); }
double TAMOMDIS_M::Sc2(double b) const{ return fSc->SMatrix2(b); }
cdouble TAMOMDIS_M::Sv(double b) const{ return fSv->SMatrix(b); }
double TAMOMDIS_M::Svr(double b) const{ return fSv->SMatrixR(b); }
double TAMOMDIS_M::Sv2(double b) const{ return fSv->SMatrix2(b); }
double TAMOMDIS_M::Svi(double b) const{ return fSv->SMatrixI(b); }

// return fMomArr[i] in MeV/c
double TAMOMDIS_M::kz(int i) const{
  if(i < 0) TAException::Error("TAMOMDIS_M", "kz: minus i");
  if(i >= kNmom) TAException::Error("TAMOMDIS_M", "kz: i=%d too large (<%d)", i, kNmom);
  return fMomArr[i]*HBARC;
} // end member function kz
